function rie = computeRie(i,elements,model)
e = model.getElement(i-1);
f = e.getFace();
B = zeros(3,8);
rie = zeros(8,1);
ng = e.getNumGaussPoint();
t = elements(i).thickness;
for j=1:ng
    xi = e.getGaussPoint(j-1);
    w = e.getGaussWeight(j-1);
    jac = f.jacXAt(xi);
    detJac = det(jac);
    B = e.computeB(f,B,xi);
    sigma = elements(i).stress(:,j);
    rie = rie + B'*sigma*(w*detJac*t);
end
if elements(i).isCrack==1
    %% compute and update renh
    elements(i) = computeRenh(e,elements(i));
    
    renh = elements(i).renh;
    rie = rie + elements(i).Kdxi*(elements(i).Kxixi\renh);
end
